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We present a new variational method, based on the matrix product operator (MPO) ansatz, for finding the 
steady state of dissipative quantum chains governed by master equations of the Lindblad form. Instead of 
requiring an accurate representation of the system evolution until the stationary state is attained, the algorithm 
directly targets the final state, thus allowing for a faster convergence when the steady state is a MPO with small 
bond dimension. Our numerical simulations for several dissipative spin models over a wide range of parameters 
illustrate the performance of the method and show that indeed the stationary state is often well described by a 
MPO of very moderate dimensions. 


Introduction .— The physics of quantum systems out of 
equilibrium poses unsolved fundamental questions, relating 
to Nature at extreme conditions and to the dynamics after 
long time evolution. Progress in this field is however hard to 
achieve, due to the lack of analytical tools to solve many such 
problems, and the limitations of existing numerical methods. 

In recent times growing attention has been directed to the 
out-of-equilibrium physics of open quantum systems, i.e. sys¬ 
tems in interaction with an environment. This interest has 
been intensified by the potential applications to the fields of 
condensed matter physics, statistical physics, and quantum in¬ 
formation processing [1-4]. In particular it has been shown 
that dissipation can be used to engineer interesting quantum 
many-body states and to perform universal quantum compu¬ 
tation [1, 5], ideas which can be explored in the context of cur¬ 
rent experimental setups based on atomic systems [6]. A par¬ 
ticularly interesting topic is that of dissipative quantum phase 
transitions (DQPT), namely transitions in the non-equilibrium 
steady state of an open system, which may arise from the com¬ 
peting effects of the Hamiltonian and the dissipative terms of 
the dynamics. An archetypical example is that of the model 
[7, 8], but DQPT have also been studied in fermionic [9-11], 
bosonic [12] and quantum spin systems [13-15], 

Finding the stationary state of a generic master equation is 
not easy, even for ID systems. Analytical treatment is limited 
to very specific problems, such as quadratic fermionic models 
[16], or systems under special conditions or approximations 
[17, 18], and most often numerical techniques are necessary. 

As in the case of pure states an exact numerical treatment 
is possible only for small systems, due to the exponentially 
growing computational cost, which may be even more severe 
in the case of mixed states. For pure states, parametrizing 
the state as a Tensor Network (TN) [19-21] has proven an ef¬ 
ficient alternative, that can successfully capture the physical 
properties of quantum many body states in countless situa¬ 
tions of interest. The best example is the tremendous success 
of the density matrix renormalization group (DMRG) [22,23], 
based on the matrix product state (MPS) ansatz, which pro¬ 
vides a quasi-exact solution for one dimensional problems. 
MPS can accurately describe ground states of gapped local 
Hamiltonians [24, 25], and methods have been defined to use 
them also in real time evolution [26-29]. In combination with 


quantum trajectories, the latter have also been applied to dis¬ 
sipative dynamics [30], The natural extension to operators, 
namely matrix product operators (MPO), can be used as an 
ansatz for mixed states [31, 32], which is known to accu¬ 
rately describe thermal equilibrium states for local Hamiltoni¬ 
ans [33, 34], Such an extension, in combination with the time 
evolution algorithms, has allowed the numerical exploration 
of steady states of spin chains and other one dimensional sys¬ 
tems under local dissipation (see e.g. [3, 35-40]). 

This method is formally similar to the search for a ground 
state using imaginary time evolution [26], in that a given ini¬ 
tial state is evolved until reaching a fixed point of the dy¬ 
namics. However, different to the imaginary time evolution 
method, where the sequence of states visited by the algorithm 
is not of physical significance, in the simulation of a mas¬ 
ter equation, the real evolution needs to be followed, so that 
errors in the intermediate state can severely affect the conver¬ 
gence of the procedure. 

A better alternative could be given by a variational method, 
which searches for the null vector of the Lindblad superoper¬ 
ator within the MPO family, in the spirit of the DMRG vari¬ 
ational search. Such a method would be potentially more ef¬ 
ficient than simulating the full evolution, specially when the 
latter traverses intermediate states with a large bond dimen¬ 
sion, but the true steady state is described by a small one, as 
is often the case [3, 36, 41], In this paper we present such 
variational method for the steady state of a master equation in 
Lindblad form. We illustrate the performance of the algorithm 
with results for several one dimensional models. Notice that 
a variational method, similar in spirit but restricted to density 
matrices containing only few-body correlations, has been re¬ 
cently proposed in [42], 

Basic concepts. — A matrix product state (MPS) for a quan¬ 
tum system of N d-dimensional components, is a state vec¬ 
tor of the form |'L) = Z{ Si } tr (A^ 1 ... A s ^ ) |si. .. sjv) [43], 
where each is a d x D x D tensor, D is a parameter of the 
representation called bond dimension, and the sum runs over 
all elements of each individual basis Si = 1,... d. By suc¬ 
cessively increasing the bond dimension, D, the MPS family 
defines a hierarchy of states covering vector space spanned by 
the tensor product of the individual bases, {|si)}. The same 
ansatz can be used to represent operators whose coefficients 
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in a tensor product basis have the structure of a matrix prod¬ 
uct, 6 = Z{si,r-i} tr (A^ iri ... A s ff rN } |si ... sjv)(t"i... rw\ 
These are called matrix product operators (MPO) [31, 32, 44]. 
The operators can be vectorized using Choi’s isomorphism, 
|si)(t*i| | SiTi), which maps any operator O to a vector 

|$(0)), so that it is possible to work in the vector space of 
operators with the usual MPS techniques. 

In order to describe physical mixed states, MPO, or in this 
case matrix product density operators (MPDO), have to sat¬ 
isfy additional conditions, namely they have to be normal¬ 
ized (trp = 1), Hermitian and positive semidefinite. While 
the first two conditions are easy to impose on the local ten¬ 
sors of the ansatz, the positivity involves the full spectmm of 
the operator, and is thus a non-local property. The ansatz can 
be modified to represent positive operators, using a local pu¬ 
rification of the state with MPS structure. In this case, each 
tensor has a structure A V = Y,k X[f ® X^ fc , where the index 
k sums over the ancillary degree of freedom and the bar indi¬ 
cates complex conjugation. Although it guarantees positivity, 
working with the purification ansatz is in general computa¬ 
tionally more costly [45] and moreover the bond dimension 
required to write the purification ansatz may be much larger 
than that of the MPO [46], so that in practice it is not always 
the most convenient choice. 

A variational search for the steady state. — We consider 
a chain of length N , with a quantum system of physical di¬ 
mension d on each site, and dynamics governed by a master 
equation of Lindblad form, ^ = C[p\. where the rhs is the 
Lindbladian superoperator, 

C[p] = -i[H,p] + £ \ (2 L aP Ll - {Li x L a ,p}) . CD 

a " 

The unitary part of the evolution is determined by the system 
Hamiltonian, H. The effect of the environment is described 
by a set of Lindblad operators, L a . 

The Lindbladian acts linearly on the vectorized p, as 

C = -i{H®t + t®H) 

+ £ i (2L a ® Z Q - ® 1 - 1 ®L T a L a ). (2) 

O' " 

The steady state is a fixed point of the evolution, - 0, 
and corresponds to a vector |$(p s )) satisfying £|$(/9 S )) = 0, 
i.e. a zero eigenvector of C. If the Hamiltonian and the in¬ 
dividual Lindblad operators have local character, the Lindbla¬ 
dian can be written as a MPO, 1 and we can search for the best 
MPS approximation to its zero eigenvector, which will give us 
a vectorized MPO approximation for the steady state. Since 
the operator (2) is not Hermitian, in order to use the standard 
variational search with MPS we consider instead the Hermi¬ 
tian product £.' C. The steady state is also a zero eigenvector 


1 Strictly speaking, it is enough that H and each L a can themselves be writ¬ 
ten as MPO, which includes short range interactions and dissipation, but 
can also be applied to approximate power-law decaying terms [44], 


of this operator, and, since t' t> 0, it corresponds to the low¬ 
est eigenvalue. If C can be written as a MPO, also the product 
can, and it is then possible to use standard MPS algorithms 
to approximate its ground state [20, 23], Notice the particular 
case of Hermitian L a is especially easy, since the (properly 
normalized) identity is a steady state, which can be exactly 
written as a MPS with bond dimension D = 1. 

The fact that we are targeting density matrices requires par¬ 
ticular attention, because not every MPS vector can represent 
a valid physical state. The normalization condition trp = 1 
translates to ($(l)|cl>(p)) = 1, where |$(1)) is the (unnor¬ 
malized) vector that corresponds to the trace map, namely the 
maximally entangled |$(1)} = Z{ Si } |si • • • sat ) <8> |si... sjv)- 
A solution which is not orthogonal to this vector can always 
be normalized to ensure the trace condition. In general it is 
more complicated to decide whether a MPS corresponds to a 
positive operator, since we do not have access to the full spec¬ 
trum. The purification ansatz can guarantee that the search 
runs over only positive operators, but at the expense of more 
costly local optimizations [45]. Hence we use simply the vec¬ 
torized MPO form and rely on the mathematical properties of 
the problem to provide a physical solution. Since the evolu¬ 
tion generated by £ is a CP map, it must have a positive fixed 
point, so that if this is non-degenerate, the algorithm should 
naturally converge towards a MPO approximation of a posi¬ 
tive (and hence Hermitian) operator 2 and is then expected to 
be almost positive, with any non-positiveness being compati¬ 
ble with the truncation error. 

In practice, we find that a suitable warmup phase (see Sup¬ 
plemental Material [48]) allows us to avoid solutions with 
vanishing trace, and improves the convergence of the algo¬ 
rithm. 3 Although positivity cannot be checked explicitly (it 
is in fact a hard problem [49]), there is a number of neces¬ 
sary criteria that any physical state needs to satisfy, such as 
physically sensible values of all single body observables. Our 
algorithm performs a set of such tests and only accepts solu¬ 
tions that pass them all, otherwise restarting the search with a 
different initial guess. The role of our consistency checks is 
to discard the least suitable guesses during the warmup phase, 
in order to prepare a suitable initial state for the variational 
search. During the later phases of the algorithm, the tests are 
used as assertions, while we rely on the convergence criteria 
(including that of the effective energy) to stop the calculation. 
After finding an acceptable solution for a given bond dimen¬ 
sion, D, we compute the desired expectation values from the 
hermitian part of the MPO, (p + p')j 2, as often done in other 


- Strictly speaking, this may fail if the algorithm gets stuck in local minima, 

or if there are degeneracies, as is known to happen also for DMRG (see 
e.g. [47]). 

3 Notice that it is also possible to use a Lagrange multiplier term, of the form 
|<f>(l))(<f > (l)|, to favor solutions with non-vanishing trace. This can be 
substracted from O C, thus increasing in one unit the bond dimension of 
the MPO. In practice, we found that the warm up phase was enough to 
obtain physical solutions, without increasing the computational cost. 
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FIG. 1. Left: Correlation (5^) for the low dimensional Dicke model, 
as a function of < 7/7 for system sizes IV = 10 - 100 (in increasingly 
darker shades). The solid line is the result of finite size extrapolation, 
linear in 1/N, as explicitly shown in the inset for several values of 
< 7 / 7 . Right: Purity of the converged steady state for the same sys¬ 
tem sizes. For large dissipation (shaded region) we show only results 
for the smallest systems, for which the simple variational search con¬ 
verges reliably. The inset shows the explicit dependence of the purity 
on the system size for several coupling values. 

algorithms to reduce numerical errors . 4 The found solution 
is normalized and used as initial guess for a larger bond di¬ 
mension, and finally convergence in D is decided when the 
targeted observables are converged to the desired precision. 

The algorithm as described here is thus formally equivalent 
to the variational ground state search for a MPO Hamiltonian 
over the MPS family, and presents the same scaling, only with 
d 2 playing the role of the physical dimension, and with an ef¬ 
fective Hamiltonian which has the squared bond dimension of 
the MPO for C. The gap of C 1 C will be determinant for the 
convergence of the algorithm. It is interesting to notice that 
this is not related to the eigenvalues of £, but to its (squared) 
singular values (see [48]). All in all we find that, for typical 
cases, the small bond dimension required to approximate the 
steady state as a MPO compensates for the additional compu¬ 
tational effort associated to £’£, provided the Lindbladian is 
not degenerate. 

Another variational approach has been recently proposed 
[42], which chooses to minimize the trace norm of £[p]. 
Since this quantity is not efficiently computable, the method 
in [42] proceeds by minimizing an upper-bound to this norm 
for restricted sets of density matrices. Our eigenvalue mini¬ 
mization is instead equivalent to finding the vector that min¬ 
imizes the Euclidean norm ||£|$(p))|| with the constraint 
|||$(p))|| = 1. Both minimizations have an exact solution 
in the physical steady state, although they are not equivalent 
when not exactly on the stationary state. Using the Euclidean 
norm of the vectorized expression is preferable in our case, 
because the trace norm requires the full diagonalization of the 
operators, impossible for the system sizes we are interested 
in, while the norm of the vectorized operators is efficiently 
calculable for the MPS ansatz. 


4 The norm of the non-Hermitian part can also be used as additional consis¬ 
tency check. 


Numerical results. — To illustrate the performance of the 
algorithm we apply it to several spin chains, where the unitary 
and dissipative dynamics show competing effects regarding 
the coherence of the steady state. 

Low dimensional Dicke model. A typical example of 
DQPT is exhibited by the Dicke model [7, 8 ], in which the 
collective interaction with a single radiation mode induces 
coherent behavior on a system of N two-level atoms. The 
regime of parameters required to observe the DQPT is chal¬ 
lenging, and the experimental observation of the phase tran¬ 
sition has only been achieved recently [50-52]. It is thus in¬ 
teresting to understand the behavior of similar models which 
may then be easier to realize experimentally. We consider a 
chain of N two-level systems, where each pair of systems 
couples coherently to a common radiation mode. This can 
be represented by a spin- 1/2 chain, governed by a single¬ 
particle Hamiltonian, H - X,= i .(/of) and Lindblad operators 
Li - 7 (<t/ + T~ +1 ), for * = 1,... TV - 1, instead of the sin¬ 
gle collective Lindblad operator of the Dicke model, so that 
this model can be considered a low-dimensional version of 
the latter. We study the nature of the steady state found by 
the algorithm at varying values of g/ 7 and increasing system 
sizes, N up to 100, which allows us to perform a finite size 
extrapolation and study single site observables and correla¬ 
tions in the thermodynamic limit. In the Dicke model, the 
superradiant phase transition (at g/ 7 2 = N) [ 8 ] is visible in 
these observables. In the low dimensional version we do not 
find evidence of such transition, although (short range) cor¬ 
relations appear, as shown in figure 1 for S ' 2 = (Zj of/TV). 
It is remarkable that for all values g/ 7 > 0.35 and system 
sizes N < 100 the steady state is converged with very small 
bond dimension, D < 30, most of them even with D < 20. 
At g = 0 there are however two dark states, namely |O) 0Ar 
and (1/s/N) £^ = 1 (-l) fc |0® (fc_ 1 h0® ( ^- fc) ). Hence, the null 
subspace of C is four-fold degenerate. This hinders the con¬ 
vergence of the algorithm at very small < 7 / 7 , as the steady 
state is no longer the unique and positive zero eigenvector, 
and the warmup strategy is not enough to guarantee a phys¬ 
ical solution, except for the smallest system sizes (N < 20 ). 
For those converged cases, we can detect the peculiarity of 
this parameter region by analyzing the purity of the solution, 
shown in figure 1. Indeed, we can find positive solutions with 
increasing purity, which can be up to 1 for g = 0 . 5 In principle 
one could complement the method with additional techniques 
to try and select the physical steady states (e.g. finding and 
then processing several orthogonal eigenstates, not necessar¬ 
ily physical) even for larger chains. Here we have nevertheless 
focused on the convergence in the most commonly occurring 
situation of a unique steady state, where the method can pro¬ 
vide the largest gain by directly targeting a MPO with small 
bond dimension. 


5 Notice that for g = 0 any mixture of both dark states will be a steady state, 
with purity in [0.5,1]. 
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FIG. 2. Left: antiferromagnetic order parameter ^/(Mf) in the 
steady state of the Ising model with local dissipation, for several 
system sizes, and fixed 7 = 1, V = 5, ft = 1.5, as a func¬ 
tion of A. The inset shows the exponential decay of correlations, 

C~J 2 (n) = {v z N/2 v z N/2+n ) - (v Z N/ 2 ){° Z N/ 2 +n)’ for N = 50 in the 

various regions of A. Right: Purity of the converged steady state for 
the same system sizes. The inset shows the local polarization, 07 , for 
the N = 10 case. At A = 0 the antiferromagnetic ordering can be ap¬ 
preciated, while for larger values of |A| the steady state approaches 
total polarization. 


Dissipative Ising chain. A complementary kind of model 
is one where the Hamiltonian dynamics induces correlations, 
for instance an Ising chain, and the dissipation is purely lo¬ 
cal. We consider a nearest neighbor Ising interaction, H = 
\ £i<i\rof 0 f + i+£i (j of - Z ^ A cr.f) + ^(crj+cr Z N ), and local 
dissipation given by Li = ^7 erf, i = 1,... IV. Such a model 
has attracted considerable attention in recent years, as it can 
be effectively realized in atomic lattice systems using Rydberg 
states [15, 53, 54]. In order to compare to existing results in 
the literature [53], we fix 7 = 1, V = 5, = 1.5. We com¬ 

pute the steady state for systems up to ./V = 50 and study the 
squared of the staggered magnetization, M z = Zj(-l) V~/7V, 
equivalent to the antiferromagnetic order parameter defined 
in [53], and the purity of the steady state as a function of A 
(Fig. 2). We find that our convergence criteria are met with 
small bond dimension D < 20. Our results show the order pa¬ 
rameter vanishing as N -*• 00 , consistent with short range cor¬ 
relations, which indeed are observed to decay exponentially. 
We observe that the purity of the steady state grows for large 
absolute values of A. This can be easily understood by going 
to the interaction picture with respect to the single body of 
terms in the Hamiltonian. This does not change the form of 
the dissipative terms, but for very large \ V - A|, the of terms 
can be neglected in the rotating wave approximation. In this 
situation the single dark state of the dissipation, the fully po¬ 
larized state 0)® ;V , is also an eigenstate of the Hamiltonian, 
and is thus a steady state. 

The method can also be applied to other models, for in¬ 
stance with coherence induced by both the Hamiltonian and 
the environment [48] . 

Conclusion. — We have presented and analyzed a varia¬ 
tional algorithm that searchs for a MPO approximation to the 
steady state of an open quantum system. The algorithm is ap¬ 
plicable to any model in which the Hamiltonian and the Lind- 
blad operators can be expressed as MPO. Instead of simulat¬ 
ing the real time evolution of the system, as done by other 


existing tensor network approaches, this method directly tar¬ 
gets the stationary state, without the need to precisely describe 
intermediate states which may need a larger bond dimension 
than the actual solution. Thus our technique can allow for a 
more efficient exploration of the steady state phase diagram. 
Our numerical results show that for a variated set of models, 
with correlations created by the unitary evolution, the dissipa¬ 
tion or both [48], the steady state is indeed well approximated 
by a MPO of very small bond dimension, D < 30 for sizes up 
to N = 100. This can be directly compared to the bond di¬ 
mensions required to describe the intermediate states in time 
evolution methods. For instance in [55] D ~ 200 was required 
for a dissipative Ising chain of length N = 40. In [3], the evo¬ 
lution required D of several hundreds [56], for a XXZ chain 
of length N = 96, when the steady state has D = 1. 

Our approach is based on the ground state optimization over 
MPS for a MPO Hamiltonian, and relies on the guaranteed 
existence of a valid, positive steady state. This basic technique 
is complemented with a warm-up phase or a suitable initial 
guess, found to be crucial in practice for convergence to a 
physical result with small bond dimension. 

When the steady state is degenerate, the simplest method 
described in this paper might have problems to find a valid 
guess for the steady state. Specially in the situation of sev¬ 
eral dark states, the null subspace of the Lindbladian contains 
infinitely many vectors which do not correspond to positive 
operators and hence do not constitute valid physical states. 6 
In principle it would be possible to complement the current 
algorithm with additional techniques, such as symmetries, in 
order to reduce the degeneracy, or to construct a candidate 
steady state from appropriate combinations of several linearly 
independent null vectors, even if non positive. 


SUPPLEMENTARY 

Basics of “ground state” searching for MPO 

In the variational search for the best MPS approximation 
to the ground state of a Hamiltonian, H , the optimization 
proceeds using an alternating least squares (ALS) strategy, in 
which the energy, (T , |fL|'I')/('I'|T'),is successively minimized 
with respect to one of the tensors in the ansatz, while the re¬ 
maining tensors are fixed [57]. In the case of the steady state 
optimization, the basic algorithm is identical, where the prod¬ 
uct C' C plays the role of the Hamiltonian. We use a MPO 
description for the superoperator C (Eq. 2 in the main text), 
which can be easily constructed as for the Hamiltonian case, 
so that also C 1 C has a MPO representation with the squared 
bond dimension of the former (see figure 3). In particular, in 
the cases described in this work, the operator C contains only 


6 Notice that this situation could also be adverse for time evolving numerical 
methods. 
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local and nearest-neighbor terms, but the construction can be 
extended to more general situations [44]. Running the ALS 
optimization requires the corresponding effective operator on 
each tensor, i.e. the contraction of ('P|£ I X|'J') except for this 
tensor, which can be computed efficiently using the same ba¬ 
sic computational routines as the ground state algorithms. 

The search for the ground state of CL can then produce a 
MPS approximation to the vector |<f>(p s )}, normalized in Eu¬ 
clidean norm, which minimizes ||£|<f>(/? s ))||. But not every 
such approximation will represent a possible physical state, 
as those need to be positive operators, with trace one. Math¬ 
ematically, there is a zero-eigenstate of t corresponding to 
a positive operator, so that the algorithm will in general be 
naturally looking for a physical solution (except in the case, 
illustrated in the main text, of several dark states). Neverthe¬ 
less we identify some additional techniques that help ensuring 
the stability and fast convergence of the algorithm. 


Initial warm-up phase. 

In general it is convenient to start the variational search with 
a small bond dimension, D, and to use the result as initial 
guess for the search with increased D until the required pre¬ 
cision is attained. In the steady state case, nevertheless, if 
the initial bond dimension is very small, it might happen that 
the algorithm at this stage converges to some vector not cor¬ 
responding to a physical state, which constitutes a bad initial 
guess for the next rounds. This may lead to very slow con¬ 
vergence, and to an artificial increase of the required bond 
dimension for the final solution. Hence, for the smallest bond 
dimension, instead of directly starting the algorithm on a ran¬ 
dom vector, we implement a first warmup phase that con¬ 
structs a more suitable initial state. This is then fed to the 
usual DMRG-like sweep, until convergence. We also observe 
that the algorithm behaves better when the bond dimension is 
increased in small steps. 

In particular, for the (close to) reflection symmetric models 
considered here, the first sweeps over the MPS tensors are not 
performed from left to right and back, as in the regular DMRG 
algorithms, but symmetrically from the outside in, or from the 
inside out. We find that this strategy works best when starting 
from the separable ansatz, D - 1, for which the algorithm is 
extremely fast, and can thus be made to converge to numer¬ 
ical precision, and be repeated if necessary at a low cost. If 
after the algorithm has converged for D = 1, if the solution is 
not compatible with a physical state, the procedure is repeated 
from a different initial state. Otherwise, the bond dimension 
is increased, and the obtained state is used as initial guess. 7 


7 When exploring the parameter space of a certain model, this technique can 
be used for some parameter values and then the converged results can be 
used as guesses for the model with slightly changed parameter. 


Physical solutions. 

Checking whether a given MPO corresponds to a physi¬ 
cal operator is in general a very hard problem. Neverthe¬ 
less, we can apply some compatibility checks to discard the 
most unphysical solutions. More concretely, for the spin-1/2 
chains studied, we demand that each of the local magnetiza¬ 
tions, tj®, <7%, erf, have expectation values within the physical 
range. s We found this to be enough to identify the most un¬ 
physical intermediate solutions, and to obtain solutions with 
non-vanishing trace, so that the steady state could be prop¬ 
erly normalized, but it is possible to define additional tests and 
perform more demanding physicality checks. Additionally, it 
might also be interesting to complement the minimization of 
C 1 L with a Lagrange multiplier term to favor vectors that are 
not orthogonal to the identity. 

Deciding convergence. 

For the ground state search, convergence of the global algo¬ 
rithm can be decided in terms of the relative variation of the 
energy from one value of the bond dimension to the next. In 
our case, the exact solution has zero eigenvalue. Instead of its 
relative variation, we check the absolute value, so we require 
that the found eigenvalue is below some threshold, and addi¬ 
tionally, demand convergence of some physical observables. 
In particular, we require that the local polarizations converge, 
by constructing the iV-component vectors, (<r“,... er^/), for 
a - x, y, z and requiring that the relative variation (in Eu¬ 
clidean norm) when increasing the bond dimension is below 
10~ 4 . With these criteria, the steady state is converged for all 
the cases presented in this work. Additionally we check that 
other observables (such as the antiferromagnetic order param¬ 
eter for the Ising model with local dissipation) and even the 
purity, which is far from being a local observable, are con¬ 
verged to a precision better than 10 -2 . 

Ising chain with coherent dissipation. 

In the main text we have illustrated the performance of 
the algorithm in two typical examples, where the Hamilto¬ 
nian and the dissipation compete to induce or destroy cor¬ 
relations. Here we complement those results with a situa¬ 
tion in which both the unitary and the dissipative dynamics 
can induce coherence. Namely, we study an Ising model, 
H = Y.i erf crf + 1 + g T,i of, with nearest neighbor Lindblad op¬ 
erators, Li = /jerf +vcr M , i < N (plus Lpf = For chains 

of up to N = 50 sites, we find convergence with remarkably 
small bond dimension, D < 20. Figure 4 illustrates the ex¬ 
pectation value of S' 2 = (Y, erf) 2 in the steady state for the 


At least within some tolerance, which can be relaxed for the smallest bond 
dimensions. 
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(a)The MPO for p is mapped to a MPS for |<E>(/9)) by Choi’s 
isomorphism. 



(b)For the Hamiltonians and Liouvillian operators considered 
in this paper, the Lindblad map, C(p) (left), is 
correspondingly mapped to a MPO acting linearly as C\Q(p)) 
(right). 



(c)The MPO sketched 
above can be used to 
compose A£|$(p)). 


FIG. 3. MPO-MPS representation of the state and the Lindblad super 
operator. 


particular case g = 0.5 and changing values of g and v, after 
extrapolating the results to TV -*• oo. It is easy to check that in 
the thermodynamic limit the identity is always a steady state 
for the symmetric case, g = v. This explains the vanishing 
(S'l) observed in the plot when u = 0.5 for all values of g. 



FIG. 4. Total 5? in the steady state for the Ising model with nearest- 
neighbor dissipation, as a function of v at fixed g = 0.5 for various 
values of g , and after finite size extrapolation from system sizes up 
to N = 50. The right plot shows the same quantity for varying g and 
fixed g = 1, for N = 40. 


£ 1 £ are given by the squared singular values of £. 

Here we illustrate how the gaps vary for the two models pre¬ 
sented in the main tex, in the case of small chains, using exact 
diagonalization. In principle it would be possible to extend 
this study and perform a finite size extrapolation, to predict 
the performance of the algorithm for increasingly larger sys¬ 
tems. Such exhaustive study escapes the scope of this paper, 
but the finite system results already allow us to appreciate the 
very different situations the algorithm can face. 

For the low dimensional Dicke model we observe that the 
gap closes, even for small systems, in the region of small 
g/ 7 . This is consistent with the parameter region on which 
we found the most difficulties. However, notice that the algo¬ 
rithm easily converged for system sizes N « 10, in spite of the 
closing gap. 

For the Ising chain with local dissipation, when we fix 
the parameters as in the text and vary only A, the gap ex¬ 
hibits much less dramatic dependence. Nevertheless, we have 
shown that in this case, easily treatable by our method, the 
steady state can present rich features. 



FIG. 5. Gap of A £ as a function of the model parameters for small 
chains, N = 3-6. For the low dimensional Dicke model (left), the 
gap closes at small values of p/ 7 . The dissipative Ising chain (right, 
for 7 = 1 , V = 5, Q. = 1.5) does not exhibit such a strong dependence 
of the gap on the free parameter, A. In both models the final values 
($(p)| A £ |4>(p)) that we have accepted are all smaller than 1 GT°. 
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Spectral gaps of £ f £ and L 

The gap of £' £ is determinant for the convergence of the 
presented algorithm. It plays the same role as the energy 
gap in conventional DMRG. For time evolution algorithms the 
most relevant quantity is the spectral gap of £, which cannot 
be related directly to the former. Instead, the eigenvalues of 
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